Load required packages
library(tidyverse)
## ── Attaching packages ─────────────────────────────────────── tidyverse 1.3.2 ──
## ✔ ggplot2 3.4.1 ✔ purrr 1.0.1
## ✔ tibble 3.1.8 ✔ dplyr 1.1.0
## ✔ tidyr 1.3.0 ✔ stringr 1.5.0
## ✔ readr 2.1.4 ✔ forcats 1.0.0
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()
library(stringi)
library(doParallel)
## Loading required package: foreach
##
## Attaching package: 'foreach'
##
## The following objects are masked from 'package:purrr':
##
## accumulate, when
##
## Loading required package: iterators
## Loading required package: parallel
library(pomp)
##
## Welcome to pomp!
##
## As of version 4.6, no user-visible pomp function has a name that
## includes a dot ('.'). Function names have been changed to replace the
## dot with an underscore ('_'). For more information, see the pomp blog:
## https://kingaa.github.io/pomp/blog.html.
##
##
## Attaching package: 'pomp'
##
## The following object is masked from 'package:purrr':
##
## map
library(panelPomp)
library(ggpubr)
Set working directory
setwd("~/Documents/GitHub/bdd/nw11_hier/final_files/")
Colour-blind friendly colour palette
cbPalette <- c("#999999", "#E69F00", "#56B4E9", "#009E73", "#F0E442", "#0072B2", "#D55E00", "#CC79A7")
Read in PNAS data
## # A tibble: 6 × 11
## day Pd RBC Ter119 CD71 mouseid paba box mouse Eryth Retic
## <int> <dbl> <dbl> <dbl> <dbl> <chr> <chr> <chr> <chr> <dbl> <dbl>
## 1 0 NA 8360000 5105135. 1.37e5 01-01 0.05 01 01 8.14e6 2.25e5
## 2 1 7934. 8290000 2745516. 1.94e5 01-01 0.05 01 01 7.70e6 5.86e5
## 3 2 19489. 7560000 2551937. 1.51e5 01-01 0.05 01 01 7.11e6 4.46e5
## 4 3 228842. 7820000 6400210. 4.98e5 01-01 0.05 01 01 7.21e6 6.08e5
## 5 4 1534425 7520000 3471975. 4.16e5 01-01 0.05 01 01 6.62e6 9.01e5
## 6 5 4528560 7600000 4748225. 5.50e5 01-01 0.05 01 01 6.72e6 8.80e5
Obtain estimates for beta and dose
## box mouse mouseid Beta dose sigmaPd sigmaRBC sigmaRetic sigmaW
## 1 01 01 01-01 6.511778 793.298642 2 0.1 0.3 1
## 2 01 02 01-02 6.511778 3.444740 2 0.1 0.3 1
## 3 01 03 01-03 6.511778 420.158396 2 0.1 0.3 1
## 4 02 01 02-01 6.090502 607.635623 2 0.1 0.3 1
## 5 02 02 02-02 6.090502 260.448444 2 0.1 0.3 1
## 6 02 03 02-03 6.090502 1.348118 2 0.1 0.3 1
## sigmaN sigmaR E_0 R_0 W_0 N_0
## 1 0.5 0.5 8e+06 3e+05 88000 8e+05
## 2 0.5 0.5 8e+06 3e+05 88000 8e+05
## 3 0.5 0.5 8e+06 3e+05 88000 8e+05
## 4 0.5 0.5 8e+06 3e+05 88000 8e+05
## 5 0.5 0.5 8e+06 3e+05 88000 8e+05
## 6 0.5 0.5 8e+06 3e+05 88000 8e+05
Create object pos with pomp object for each mouse
foreach (m = iter(theta,"row"),.inorder=TRUE,.combine=c) %dopar% {
flow %>%
filter(mouseid==m$mouseid) %>%
select(day,Pd,RBC,Retic) %>%
mutate(Retic=if_else(day %in% c(0,14),NA_real_,Retic)) %>%
pomp(
params=select(m,-mouseid,-box,-mouse) %>% unlist(),
times="day",
t0=0,
rmeasure=Csnippet("
Retic = rlnorm(log(1+R),sigmaRetic)-1;
RBC = rlnorm(log(1+E+R),sigmaRBC)-1;
Pd = rlnorm(log(1+K),sigmaPd)-1;"),
dmeasure=Csnippet("
double l1, l2, l3;
l1 = (R_FINITE(Retic)) ? dlnorm(1+Retic,log(1+R),sigmaRetic,1) : 0;
l2 = (R_FINITE(RBC)) ? dlnorm(1+RBC,log(1+E+R),sigmaRBC,1) : 0;
l3 = (R_FINITE(Pd) && Pd>0) ? dlnorm(1+Pd,log(1+K),sigmaPd,1) : 0;
lik = (give_log) ? l1+l2+l3 : exp(l1+l2+l3);"),
rprocess=discrete_time(
step.fun=Csnippet("
double Mold = M;
M = Beta*K*exp(-(W+N)/(R+E));
E = (R+E)*exp(-(Mold+N)/(R+E));
N = rlnorm(log(N),sigmaN);
W = rlnorm(log(W),sigmaW);
R = rlnorm(log(R),sigmaR);
K = (R+E>0) ? (R+E)*(1-exp(-M/(R+E))): 0;
"),
delta.t=1
),
partrans=parameter_trans(
log=c("sigmaW","sigmaR","sigmaN",
"sigmaPd","sigmaRBC","sigmaRetic",
"N_0","W_0","E_0","R_0")
),
rinit=Csnippet("
E = E_0;
R = R_0;
N = N_0;
W = W_0;
M = 0;
K = dose;"),
statenames=c("E","R","W","N","M","K"),
paramnames=c(
"Beta","dose",
"sigmaPd","sigmaRBC","sigmaRetic",
"sigmaW","sigmaN","sigmaR",
"E_0","R_0","W_0","N_0"
)
)
} %>%
set_names(theta$mouseid) -> pos
Obtain MLE for sigmas, initial values, betas and dose
## loglik sigmaW sigmaN sigmaPd sigmaR sigmaRBC sigmaRetic E_0
## 1 -11217.79 1.043139 0.624558 0.5966612 0.392552 0.07942614 0.1514119 7722711
## N_0 R_0 W_0
## 1 410608.7 617657.5 176481.2
## # A tibble: 15 × 3
## mouseid Beta dose
## <chr> <dbl> <dbl>
## 1 01-01 6.51 793.
## 2 01-02 6.51 3.44
## 3 01-03 6.51 420.
## 4 02-01 6.09 608.
## 5 02-02 6.09 260.
## 6 02-03 6.09 1.35
## 7 03-01 5.76 591.
## 8 03-02 5.76 382.
## 9 03-03 5.76 317.
## 10 04-01 4.08 1476.
## 11 04-02 4.08 716.
## 12 04-03 4.08 1139.
## 13 05-01 0 0
## 14 05-02 0 0
## 15 05-03 0 0
Load in PNAS trajectories
sm1name <- "m5sm1.rds"
sm1 <- readRDS(sm1name)
Create dataframe with weighted trajectories (grouped by mouseid
first, then by box)
sm1 |>
as_tibble() |>
filter(mouseid!="01-02",mouseid!="02-03") |> #remove underdosed mice
pivot_wider(names_from=variable,values_from=value) |>
left_join(bdf,by=c("mouseid")) |>
separate_wider_delim(cols="mouseid",delim="-",names=c("box","mouse"),cols_remove=FALSE) |>
group_by(mouseid) |>
mutate(
SM=exp(-M/(R+E)),
SN=exp(-N/(R+E)),
SW=exp(-W/(R+E)),
Qps=(1-SM)*SW*SN,
Qpn=N/(N+W)*(1-SM)*(1-SW*SN),
Qpw=W/(N+W)*(1-SM)*(1-SW*SN),
Qun=SM*(1-SN),
Qus=SM*SN,
lambda_r=Beta*R/M*Qps,
lambda_e=Beta*E/M*Qps,
lambda_n=Beta*(R+E)/M*Qpn,
lambda_w=Beta*(R+E)/M*Qpw,
lambda_u=Beta-lambda_r-lambda_e-lambda_n-lambda_w,
lambda_u_w_ratio=lambda_u/lambda_w,
lambda_u_n_ratio=lambda_u/lambda_n,
lambda_u_wn_ratio=lambda_u/(lambda_w+lambda_n),
loss=(R+E)*(1-Qus),
rbc=E+R,
varN=N/rbc,
perRetic=R/rbc,
lik=exp(loglik-max(loglik))
) |>
ungroup() |>
select(-loglik,-Beta,-dose) |>
gather(variable,value,-rep,-time,-mouse,-box,-lik,-mouseid) |>
filter(is.finite(value)) %>%
group_by(box,time,variable) |>
dplyr::reframe(
value=pomp::wquant(x=value,probs=c(0.05,0.5,0.95),weights=lik),
name=c("lo","med","hi")
) |>
ungroup() |>
pivot_wider() -> group_traj
group_traj$pABA <- factor(group_traj$box,levels=c("05","04","03","02","01"),
labels=c("Control","0%","0.05%","0.5%","5%"))
Remove estimates for W for control mice
group_traj <- group_traj |> dplyr::slice(-which(group_traj$box=="05"&group_traj$variable=="W"))
group_traj |>
filter(time<=20,variable%in%c("rbc","E","R","N","W","K")) |>
mutate(variable=case_match(variable,
"K"~"Parasites",
"rbc"~"RBC",
"R"~"Reticulocytes",
"E"~"Erythrocytes",
"N"~"Indiscriminate killing",
"W"~"Targeted killing")) |>
ggplot()+
geom_line(aes(x=time,y=med,col=pABA))+
geom_ribbon(aes(x=time,ymin=lo,ymax=hi,fill=pABA),alpha=0.2)+
scale_y_log10()+
scale_colour_manual(values=cbPalette)+
scale_fill_manual(values=cbPalette)+
xlab("Days")+ylab("Density per microlitre")+
facet_wrap(variable~.,scales="free_y")+
theme_bw()+
theme(
strip.background = element_blank()
)
## Warning: Transformation introduced infinite values in continuous y-axis
## Transformation introduced infinite values in continuous y-axis
## Transformation introduced infinite values in continuous y-axis
Plot group-level trajectories for reticulocytes with day of peak
parasitaemia
Plot group-level trajectories for parasite density
## Warning: Transformation introduced infinite values in continuous y-axis
## Transformation introduced infinite values in continuous y-axis
## Transformation introduced infinite values in continuous y-axis
Plot group-level trajectories for percent reticulocytes
Create swirly data frame
## # A tibble: 6 × 29
## box time pABA E K M N Qpn Qps Qpw Qun
## <chr> <int> <fct> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 01 0 5% 7722711. 4.20e2 0 4.11e5 0 0 0 0.0480
## 2 01 1 5% 7939703. 2.55e3 2.55e3 4.44e5 1.52e-5 2.78e-4 1.02e-5 0.0503
## 3 01 2 5% 8165549. 1.56e4 1.56e4 4.74e5 9.34e-5 1.63e-3 1.04e-4 0.0532
## 4 01 3 5% 8166517. 9.22e4 9.27e4 4.66e5 5.29e-4 9.04e-3 1.15e-3 0.0523
## 5 01 4 5% 8085821. 5.03e5 5.18e5 4.40e5 2.45e-3 4.01e-2 1.54e-2 0.0470
## 6 01 5 5% 7717184. 2.01e6 2.29e6 4.26e5 7.74e-3 7.99e-2 1.48e-1 0.0370
## # … with 18 more variables: Qus <dbl>, R <dbl>, SM <dbl>, SN <dbl>, SW <dbl>,
## # W <dbl>, loss <dbl>, perRetic <dbl>, rbc <dbl>, varN <dbl>, lambda_e <dbl>,
## # lambda_n <dbl>, lambda_r <dbl>, lambda_u <dbl>, lambda_u_n_ratio <dbl>,
## # lambda_u_w_ratio <dbl>, lambda_u_wn_ratio <dbl>, lambda_w <dbl>
Parasite density versus targeted killing
## Warning: Transformation introduced infinite values in continuous x-axis
## Transformation introduced infinite values in continuous x-axis
## Warning: Removed 7 rows containing missing values (`geom_path()`).
## Warning: Removed 2 rows containing missing values (`geom_text()`).
## Warning: Transformation introduced infinite values in continuous x-axis
## Transformation introduced infinite values in continuous x-axis
## Warning: Removed 3 rows containing missing values (`geom_path()`).
## Warning: Removed 3 rows containing missing values (`geom_text()`).
## Warning: Transformation introduced infinite values in continuous x-axis
## Transformation introduced infinite values in continuous x-axis
## Warning: Removed 14 rows containing missing values (`geom_path()`).
## Warning: Removed 2 rows containing missing values (`geom_text()`).
## Warning: Transformation introduced infinite values in continuous x-axis
## Transformation introduced infinite values in continuous x-axis
## Warning: Removed 21 rows containing missing values (`geom_path()`).
## Warning: Removed 5 rows containing missing values (`geom_text()`).
## Warning: Transformation introduced infinite values in continuous x-axis
## Transformation introduced infinite values in continuous x-axis
## Warning: Removed 21 rows containing missing values (`geom_path()`).
## Warning: Removed 5 rows containing missing values (`geom_text()`).
## Warning: Transformation introduced infinite values in continuous x-axis
## Transformation introduced infinite values in continuous x-axis
## Transformation introduced infinite values in continuous x-axis
## Transformation introduced infinite values in continuous x-axis
## Transformation introduced infinite values in continuous x-axis
## Warning: Removed 3 rows containing missing values (`geom_errorbarh()`).
## Warning: Removed 3 rows containing missing values (`geom_path()`).
## Warning: Removed 3 rows containing missing values (`geom_text()`).
Indiscriminate killing versus reticulocyte supply
Targeted killing versus reticulocyte supply
## Warning: Removed 7 rows containing missing values (`geom_path()`).
## Warning: Removed 2 rows containing missing values (`geom_text()`).
## Warning: Removed 14 rows containing missing values (`geom_path()`).
## Warning: Removed 2 rows containing missing values (`geom_text()`).
## Warning: Removed 21 rows containing missing values (`geom_path()`).
## Warning: Removed 5 rows containing missing values (`geom_text()`).
## Warning: Removed 21 rows containing missing values (`geom_path()`).
## Warning: Removed 5 rows containing missing values (`geom_text()`).
Reticulocyte supply versus erythrocyte density
Plot N versus R with bi-directional error bars
Ratio of \(lambda_u\) to \(lambda_w\) over time
## Warning in self$trans$transform(x): NaNs produced
## Warning: Transformation introduced infinite values in continuous y-axis
Ratio of \(lambda_u\) to \(lambda_w\)+\(lambda_n\) over time
## Warning in self$trans$transform(x): NaNs produced
## Warning: Transformation introduced infinite values in continuous y-axis
Read in data of weighted individual trajectories
## Rows: 606492 Columns: 11
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr (2): box, mouse
## dbl (9): rep, time, E, R, W, N, K, lik, rep2
##
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
Plot individual weighted trajectories with group-level
trajectories
## Warning: Transformation introduced infinite values in continuous x-axis
## Warning: Transformation introduced infinite values in continuous y-axis
## Warning: Transformation introduced infinite values in continuous x-axis
## Warning: Removed 26 rows containing missing values (`geom_path()`).
## Warning: Transformation introduced infinite values in continuous y-axis
## Removed 26 rows containing missing values (`geom_path()`).
